Analysis of GLDS-38 from NASA GeneLab

This R markdown file was auto-generated by the iDEP website Using iDEP 0.91, originally by Steven Xijin.Ge@sdstate.edu

Ge SX, Son EW, Yao R: iDEP: an integrated web application for differential expression and pathway analysis of RNA-Seq data. BMC Bioinformatics 2018, 19(1):534. PMID:30567491

1. Read data

First we set up the working directory to where the files are saved.

input_biclustMethod_ <- "BCCC()"

R packages used

library(RSQLite, verbose = FALSE) # for database connection
library(gplots, verbose = FALSE) # for hierarchical clustering
library(ggplot2, verbose = FALSE) # graphics
library(e1071, verbose = FALSE) # computing kurtosis
#library(DT, verbose = FALSE) # for renderDataTable
library(plotly, verbose = FALSE) # for interactive heatmap
library(reshape2, verbose = FALSE) # for melt correlation matrix in heatmap

# From Data Read Function
library(edgeR, verbose = FALSE) # count data D.E.
library(DESeq2, verbose = FALSE) # count data analysis, DEG.DESeq2

# TSNE Plot, tSNEgenePlot
library(Rtsne, verbose = FALSE)

# PGSA Pathway PGSEA Pathway, PGSEAplot
library(PGSEA, verbose = FALSE)

# DEG.limma
library(limma, verbose = FALSE) # Differential expression
library(statmod, verbose = FALSE)

# enrichment plot
library(dendextend) # customizing tree

# enrich.net2, moduleNetwork
library(igraph)

# Stringdb_geneList, StringDB_GO_enrichmentData, stringDB_network1
# StringDB_network_link
library(STRINGdb, verbose = FALSE)

# gagePathwayData
library(gage, verbose = FALSE) # pathway analysis

# fgseaPathwayData
library(fgsea, verbose = FALSE) # fast GSEA

# ReactomePAPathwayData
library(ReactomePA, verbose = FALSE) # pathway analysis

# KeggImage
library(pathview)

# genomePlot, genomePlotDataPre
library(PREDA, verbose = FALSE) # showing expression on genome
library(PREDAsampledata, verbose = FALSE)
library(hgu133plus2.db, verbose = FALSE)

# biclustering
library(biclust, verbose = FALSE)

library(knitr) #  install if needed. for showing tables with kable
library(kableExtra)

if (input_biclustMethod_ == "BCQU()") {
  library(QUBIC, verbose = FALSE)
} # have trouble installing on Linux
if (input_biclustMethod_ == "BCUnibic()") {
  library(runibic, verbose = FALSE)
} # Test biclustMethod dependant qubic runibic

# wgcna
library(WGCNA)
library(flashClust, verbose = FALSE)
source("iDEP_core_functions_only.R")
# Each row of this matrix represents a color scheme;
mycolors_ <- sort(rainbow(20))[c(1, 20, 10, 11, 2, 19, 3, 12, 4, 13, 5, 14, 6, 15, 7, 16, 8, 17, 9, 18)]

hmcols_ <- colorRampPalette(colors = c('#4575B4', '#91BFDB', '#E0F3F8', '#FFFFBF', '#FEE090', '#FC8D59', '#D73027'))(75)

heatColors_ <- rbind(
  greenred(75), 
  bluered(75), 
  colorpanel(75, "green", "black", "magenta"), 
  colorpanel(75, "blue", "yellow", "red"), 
  hmcols_
)

rownames(heatColors_) <- c("Green-Black-Red", "Blue-White-Red", "Green-Black-Magenta", "Blue-Yellow-Red", "Blue-white-brown")

We are using the downloaded gene expression file where gene IDs has been converted to Ensembl gene IDs. This is because the ID conversion database is too large to download. You can use your original file if your file uses Ensembl ID, or you do not want to use the pathway files available in iDEP (or it is not available).

inputFolderFiles <- list.files(params$input_folder, full.names = TRUE)

inputFile_ <- inputFolderFiles[stringr::str_detect(tolower(inputFolderFiles), "expression.csv")]
sampleInfoFile_ <- inputFolderFiles[stringr::str_detect(tolower(inputFolderFiles), "sampleinfo.csv")]
gldsMetadataFile_ <- inputFolderFiles[stringr::str_detect(tolower(inputFolderFiles), "metadata.csv")]


geneInfoFile_ <- params$geneInfoFile
geneSetFile_ <- params$geneSetFile # pathway database in SQL; can be GMT format

STRING10_speciesFile_ <- "https://raw.githubusercontent.com/iDEP-SDSU/idep/master/shinyapps/idep/STRING10_species.csv"
readMetadata.out_ <- readMetadata(inFile = gldsMetadataFile_) #gldsMetadataFile_)

kable(readMetadata.out_) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
FLT_Rep1 FLT_Rep2 FLT_Rep3 GC_Rep1 GC_Rep2 GC_Rep3 LN2_Rep1 LN2_Rep2 LN2_Rep3 RNAlat_Rep1 RNAlat_Rep2 RNAlat_Rep3
Sample.LongId Atha.WT.Col.0.sl.FLT.Rep1.G1S1.RNAseq.RNAseq Atha.WT.Col.0.sl.FLT.Rep2.G1S2.RNAseq.RNAseq Atha.WT.Col.0.sl.FLT.Rep3.G1S3.RNAseq.RNAseq Atha.WT.Col.0.sl.GC.Rep1.G2S1.RNAseq.RNAseq Atha.WT.Col.0.sl.GC.Rep2.G2S2.RNAseq.RNAseq Atha.WT.Col.0.sl.GC.Rep3.G2S3.RNAseq.RNAseq Atha.WT.Col.0.sl.LN2.Rep1.n2.1.RNAseq.RNAseq Atha.WT.Col.0.sl.LN2.Rep2.n2.2.RNAseq.RNAseq Atha.WT.Col.0.sl.LN2.Rep3.n2.3.RNAseq.RNAseq Atha.WT.Col.0.sl.RNAlat.Rep1.rl1.RNAseq.RNAseq Atha.WT.Col.0.sl.RNAlat.Rep2.rl2.RNAseq.RNAseq Atha.WT.Col.0.sl.RNAlat.Rep3.rl3.RNAseq.RNAseq
Sample.Id Atha.WT.Col.0.sl.FLT.Rep1.G1S1 Atha.WT.Col.0.sl.FLT.Rep2.G1S2 Atha.WT.Col.0.sl.FLT.Rep3.G1S3 Atha.WT.Col.0.sl.GC.Rep1.G2S1 Atha.WT.Col.0.sl.GC.Rep2.G2S2 Atha.WT.Col.0.sl.GC.Rep3.G2S3 Atha.WT.Col.0.sl.LN2.Rep1.n2.1 Atha.WT.Col.0.sl.LN2.Rep2.n2.2 Atha.WT.Col.0.sl.LN2.Rep3.n2.3 Atha.WT.Col.0.sl.RNAlat.Rep1.rl1 Atha.WT.Col.0.sl.RNAlat.Rep2.rl2 Atha.WT.Col.0.sl.RNAlat.Rep3.rl3
Sample.Name Atha_WT-Col-0_sl_FLT_Rep1_G1S1 Atha_WT-Col-0_sl_FLT_Rep2_G1S2 Atha_WT-Col-0_sl_FLT_Rep3_G1S3 Atha_WT-Col-0_sl_GC_Rep1_G2S1 Atha_WT-Col-0_sl_GC_Rep2_G2S2 Atha_WT-Col-0_sl_GC_Rep3_G2S3 Atha_WT-Col-0_sl_LN2_Rep1_n2-1 Atha_WT-Col-0_sl_LN2_Rep2_n2-2 Atha_WT-Col-0_sl_LN2_Rep3_n2-3 Atha_WT-Col-0_sl_RNAlat_Rep1_rl1 Atha_WT-Col-0_sl_RNAlat_Rep2_rl2 Atha_WT-Col-0_sl_RNAlat_Rep3_rl3
GLDS 38 38 38 38 38 38 38 38 38 38 38 38
Accession GLDS-38 GLDS-38 GLDS-38 GLDS-38 GLDS-38 GLDS-38 GLDS-38 GLDS-38 GLDS-38 GLDS-38 GLDS-38 GLDS-38
Hardware BRIC BRIC BRIC BRIC BRIC BRIC BRIC BRIC BRIC BRIC BRIC BRIC
Tissue Etiolated seedling Etiolated seedling Etiolated seedling Etiolated seedling Etiolated seedling Etiolated seedling Etiolated seedling Etiolated seedling Etiolated seedling Etiolated seedling Etiolated seedling Etiolated seedling
Age 8 days 8 days 8 days 8 days 8 days 8 days 8 days 8 days 8 days 8 days 8 days 8 days
Organism Arabidopsis thaliana Arabidopsis thaliana Arabidopsis thaliana Arabidopsis thaliana Arabidopsis thaliana Arabidopsis thaliana Arabidopsis thaliana Arabidopsis thaliana Arabidopsis thaliana Arabidopsis thaliana Arabidopsis thaliana Arabidopsis thaliana
Ecotype Col-0 Col-0 Col-0 Col-0 Col-0 Col-0 Col-0 Col-0 Col-0 Col-0 Col-0 Col-0
Genotype WT WT WT WT WT WT WT WT WT WT WT WT
Variety Col-0 WT Col-0 WT Col-0 WT Col-0 WT Col-0 WT Col-0 WT Col-0 WT Col-0 WT Col-0 WT Col-0 WT Col-0 WT Col-0 WT
Radiation Cosmic radiation Cosmic radiation Cosmic radiation Background Earth Background Earth Background Earth Background Earth Background Earth Background Earth Background Earth Background Earth Background Earth
Gravity Microgravity Microgravity Microgravity Terrestrial Terrestrial Terrestrial Terrestrial Terrestrial Terrestrial Terrestrial Terrestrial Terrestrial
Developmental Etiolated 8 day old seedlings Etiolated 8 day old seedlings Etiolated 8 day old seedlings Etiolated 8 day old seedlings Etiolated 8 day old seedlings Etiolated 8 day old seedlings Etiolated 8 day old seedlings Etiolated 8 day old seedlings Etiolated 8 day old seedlings Etiolated 8 day old seedlings Etiolated 8 day old seedlings Etiolated 8 day old seedlings
Time.series.or.Concentration.gradient Single time point Single time point Single time point Single time point Single time point Single time point Single time point Single time point Single time point Single time point Single time point Single time point
Light Dark Dark Dark Dark Dark Dark Dark Dark Dark Dark Dark Dark
Assay..RNAseq. RNAseq Transcription and Proteomic Profiling RNAseq Transcription and Proteomic Profiling RNAseq Transcription and Proteomic Profiling RNAseq Transcription and Proteomic Profiling RNAseq Transcription and Proteomic Profiling RNAseq Transcription and Proteomic Profiling RNAseq Transcription and Proteomic Profiling RNAseq Transcription and Proteomic Profiling RNAseq Transcription and Proteomic Profiling RNAseq Transcription and Proteomic Profiling RNAseq Transcription and Proteomic Profiling RNAseq Transcription and Proteomic Profiling
Temperature 22-24 22-24 22-24 22-24 22-24 22-24 22-24 22-24 22-24 22-24 22-24 22-24
Treatment.type Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity Proteomics and Transcriptomics analysis of Arabidopsis Seedlings in Microgravity
Treatment.intensity X X X X X X X X X X X X
Treament.timing X X X X X X X X X X X X
Preservation.Method. RNAlater RNAlater RNAlater RNAlater RNAlater RNAlater Liquid Nitrogen Liquid Nitrogen Liquid Nitrogen RNAlater RNAlater RNAlater 
readData.out_ <- readData(inFile = inputFile_, 
                          input_missingValue = "geneMedian", 
                          input_dataFileFormat = 1, 
                          input_minCounts = 0.5, 
                          input_NminSamples = 1, 
                          input_countsLogStart = 4, 
                          input_CountsTransform = 1)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
kable(head(readData.out_$data)) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
FLT_Rep1 FLT_Rep2 FLT_Rep3 GC_Rep1 GC_Rep2 GC_Rep3 LN2_Rep1 LN2_Rep2 LN2_Rep3 RNAlat_Rep1 RNAlat_Rep2 RNAlat_Rep3
ATCG00490 19.02600 20.20307 20.47902 19.35267 20.26267 19.81133 17.85944 18.14881 18.28268 19.90399 19.24373 18.83189
AT2G41310 18.59846 19.97237 19.85194 20.18207 18.98753 19.11273 18.29485 18.20631 18.20418 18.69872 18.47971 18.18740
ATCG00020 17.76784 18.79161 19.31738 17.92858 19.00300 18.51242 17.11173 17.35404 17.37928 18.97660 18.68141 18.13132
AT3G21720 17.38043 17.29131 17.40174 18.45645 16.52706 17.41765 16.35671 15.58054 15.69393 14.72961 14.46187 14.90556
AT2G07671 16.80806 17.55228 17.74476 16.89524 17.38578 17.32608 15.95185 15.89930 15.66004 16.40474 16.05359 15.88383
ATCG00280 16.09754 17.02548 17.38910 17.07244 17.48614 17.20404 15.47935 15.34283 15.45231 17.24330 16.73648 16.46234 
readSampleInfo.out_ <- readSampleInfo(inFile = sampleInfoFile_, 
                                      readData.out = readData.out_)
kable(readSampleInfo.out_) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
Gravity Preservation.Method.
FLT_Rep1 Microgravity RNAlater
FLT_Rep2 Microgravity RNAlater
FLT_Rep3 Microgravity RNAlater
GC_Rep1 Terrestrial RNAlater
GC_Rep2 Terrestrial RNAlater
GC_Rep3 Terrestrial RNAlater
LN2_Rep1 Terrestrial Liquid Nitrogen
LN2_Rep2 Terrestrial Liquid Nitrogen
LN2_Rep3 Terrestrial Liquid Nitrogen
RNAlat_Rep1 Terrestrial RNAlater
RNAlat_Rep2 Terrestrial RNAlater
RNAlat_Rep3 Terrestrial RNAlater 
input_noIDConversion_ <- TRUE

allGeneInfo.out_ <- geneInfo(fileName = geneInfoFile_)
converted.out_ <- NULL
convertedData.out_ <- convertedData(converted.out = NULL, 
                                    readData.out = readData.out_, 
                                    input_noIDConversion = TRUE)

nGenesFilter(readData.out = readData.out_, 
             converted.out = NULL, 
             convertedData.out = convertedData.out_, 
             input_noIDConversion = TRUE)
## [1] "16156 genes in 12 samples. 16117  genes passed filter.\n Original gene IDs used."
convertedCounts.out_ <- convertedCounts(readData.out = readData.out_, converted.out = NULL) # converted counts, just for compatibility

2. Pre-process

# Read counts per library
parDefault_ <- par()
par(mar = c(12, 4, 2, 2))
# barplot of total read counts
rawCounts <- readData.out_$rawCounts
groups_ <- as.factor(detectGroups(colnames(rawCounts)))
if (nlevels(groups_) <= 1 | nlevels(groups_) > 20) {
  col1_ <- "green"
} else {
  col1_ <- rainbow(nlevels(groups_))[groups_]
}

barplot(colSums(readData.out_$rawCounts) / 1e6,
  col = col1_, las = 3, main = "Total read counts (millions)"
)

readCountsBias(readData.out = readData.out_, readSampleInfo.out = readSampleInfo.out_) # detecting bias in sequencing depth
## [1] 0.02579404
## [1] 0.1547212
## [1] 0.05122098
## [1] "Warning! Sequencing depth bias detected. Total read counts are significantly different among sample groups (p= 2.58e-02 ) based on ANOVA."
# Box plot
boxplot(
  x = readData.out_$data,
  las = 2, col = col1_,
  ylab = "Transformed expression levels",
  main = "Distribution of transformed data"
)

# Density plot
par(parDefault_)
## Warning in par(parDefault_): graphical parameter "cin" cannot be set
## Warning in par(parDefault_): graphical parameter "cra" cannot be set
## Warning in par(parDefault_): graphical parameter "csi" cannot be set
## Warning in par(parDefault_): graphical parameter "cxy" cannot be set
## Warning in par(parDefault_): graphical parameter "din" cannot be set
## Warning in par(parDefault_): graphical parameter "page" cannot be set
densityPlot(readData.out = readData.out_, 
            mycolors = mycolors_)

# Scatter plot of the first two samples
plot(
  x = readData.out_$data[, 1:2],
  xlab = colnames(readData.out_$data)[1], 
  ylab = colnames(readData.out_$data)[2],
  main = "Scatter plot of first two samples"
)

#### plot gene or gene family
genePlot(allGeneInfo.out = allGeneInfo.out_, 
         convertedData.out = convertedData.out_, 
         input_selectOrg = "BestMatch", 
         input_geneSearch = "HOXA")
## NULL
geneBarPlotError(convertedData.out = convertedData.out_, 
                 allGeneInfo.out = allGeneInfo.out_, 
                 input_selectOrg = 'BestMatch', 
                 input_geneSearch = "HOXA", 
                 input_useSD = "FALSE") # Use standard deviation instead of standard error in error bar?
## NULL

3. Heatmap

# hierarchical clustering tree
x <- readData.out_$data
maxGene <- apply(x, 1, max)
# remove bottom 25% lowly expressed genes, which inflate the PPC
x <- x[which(maxGene > quantile(maxGene)[1]), ]
plot(as.dendrogram(hclust2(dist2(t(x)))), ylab = "1 - Pearson C.C.", type = "rectangle")

# Correlation matrix
#input_labelPCC_ <- TRUE # Show correlation coefficient?
correlationMatrix(readData.out = readData.out_, input_labelPCC = TRUE)

png("heatmap.png", width = 10, height = 15, units = "in", res = 300)
staticHeatmap(readData.out = readData.out_, 
              readSampleInfo.out = readSampleInfo.out_, 
              heatColors = heatColors_, 
              input_nGenes = 1000, 
              input_geneCentering = TRUE, 
              input_sampleCentering = FALSE, 
              input_geneNormalize = FALSE, 
              input_sampleNormalize = FALSE, 
              input_noSampleClustering = FALSE, 
              input_heatmapCutoff = 4, 
              input_distFunctions = 1, 
              input_hclustFunctions = 1, 
              input_heatColors1 = 1, 
              input_selectFactorsHeatmap = 'Gravity')
dev.off()
## quartz_off_screen 
##                 2

[heatmap] (heatmap.png)

heatmapPlotly(convertedData.out = convertedData.out_, 
              heatColors = heatColors_, 
              allGeneInfo.out = allGeneInfo.out_, 
              input_geneCentering = TRUE, 
              input_sampleCentering = FALSE, 
              input_geneNormalize = FALSE, 
              input_sampleNormalize = FALSE, 
              input_heatColors1 = 1)# interactive heatmap using Plotly

4. K-means clustering

distributionSD(convertedData.out = convertedData.out_, 
               input_nGenesKNN = 2000) # Distribution of standard deviations

KmeansNclusters(convertedData.out = convertedData.out_, 
                input_nGenesKNN = 2000) # Number of clusters

Kmeans.out_ <- Kmeans(convertedData.out = convertedData.out_,
                      maxGeneClustering = 12000, 
                      input_nGenesKNN = 2000, 
                      input_nClusters = 4, 
                      input_kmeansNormalization = "geneMean", 
                      input_KmeansReRun = 0) # Running K-means

KmeansHeatmap(Kmeans.out = Kmeans.out_, 
              .mycolors = mycolors_, 
              .heatColors = heatColors_, 
              .input_heatColors1 = 1) # Heatmap for k-Means

# Read gene sets for enrichment analysis
GeneSets.out_ <- readGeneSets(
  fileName = geneSetFile_,
  convertedData = convertedData.out_, 
  GO = "GOBP", 
  selectOrg = "NEW",
  myrange = c(15, 2000)
)

# Alternatively, users can use their own GMT files by
# GeneSets.out_ <- readGMTRobust('somefile.GMT')
results <- KmeansGO(Kmeans.out = Kmeans.out_, 
                    input_nClusters = 4, 
                    GeneSets.out = GeneSets.out_) # Enrichment analysis for k-Means clusters

results$adj.Pval <- format(results$adj.Pval, digits = 3)

kable(results, row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
Cluster adj.Pval Genes Pathways
A 7.99e-122 219 Organonitrogen compound biosynthetic process
2.84e-118 155 Translation
4.26e-118 155 Peptide biosynthetic process
1.81e-116 159 Amide biosynthetic process
3.07e-113 156 Peptide metabolic process
1.22e-106 160 Cellular amide metabolic process
5.05e-51 157 Response to abiotic stimulus
4.90e-42 57 Photosynthesis
5.59e-35 50 Response to cytokinin
7.30e-34 63 Generation of precursor metabolites and energy
B 9.12e-48 81 Cell wall organization or biogenesis
9.12e-48 73 Cell wall organization
2.49e-47 74 External encapsulating structure organization
1.20e-40 128 Response to abiotic stimulus
2.58e-36 74 Drug metabolic process
8.64e-36 111 Small molecule metabolic process
2.30e-34 82 Carbohydrate metabolic process
2.27e-32 77 Response to inorganic substance
6.60e-28 42 Plant-type cell wall organization or biogenesis
1.13e-27 51 Polysaccharide metabolic process
C 3.98e-50 79 Response to external stimulus
1.34e-42 63 Response to external biotic stimulus
1.34e-42 63 Response to other organism
2.33e-42 63 Response to biotic stimulus
4.56e-40 42 Secondary metabolic process
1.25e-38 63 Defense response
2.03e-38 67 Multi-organism process
2.35e-35 78 Response to abiotic stimulus
4.27e-34 26 Indole-containing compound metabolic process
1.99e-33 49 Defense response to other organism
D 1.20e-50 151 Response to abiotic stimulus
1.95e-47 127 Response to oxygen-containing compound
2.73e-46 98 Response to inorganic substance
1.10e-39 101 Response to acid chemical
2.11e-39 113 Cellular response to chemical stimulus
6.77e-37 126 Response to organic substance
1.17e-31 107 Response to hormone
4.57e-31 107 Response to endogenous stimulus
3.63e-27 75 Response to external biotic stimulus
3.63e-27 75 Response to other organism 
tSNEgenePlot(Kmeans.out_, 
             input_seedTSNE = 0, 
             input_colorGenes = TRUE, # Color genes in t-SNE plot?
             mycolors = mycolors_) # Plot genes using t-SNE

5. PCA and beyond

PCAplot(convertedData.out = convertedData.out_, 
        readSampleInfo.out = readSampleInfo.out_, 
        input_selectFactors = "Gravity", 
        input_selectFactors2 = "Preservation.Method.")

MDSplot(convertedData.out = convertedData.out_, 
        readSampleInfo.out = readSampleInfo.out_, 
        input_selectFactors = "Gravity", 
        input_selectFactors2 = "Preservation.Method.")

tSNEplot(convertedData.out = convertedData.out_, 
         readSampleInfo.out = readSampleInfo.out_, 
         input_selectFactors = "Gravity", 
         input_selectFactors2 = "Preservation.Method.", 
         input_tsneSeed2 = 0)

# Read gene sets for pathway analysis using PGSEA on principal components
GeneSets.out_ <- readGeneSets(
  fileName = geneSetFile_,
  convertedData = convertedData.out_, 
  GO = "GOBP", 
  selectOrg = "NEW",
  myrange = c(15, 2000)
)

PCApathway(convertedData.out = convertedData.out_, 
           GeneSets.out = GeneSets.out_) # Run PGSEA analysis

cat(
  PCA2factor(readData.out = readData.out_, 
             readSampleInfo.out = readSampleInfo.out_)
) # The correlation between PCs with factors
## 
##  Correlation between Principal Components (PCs) with factors
## PC1 is correlated with Gravity (p=3.91e-02).

6. DEG1

limma.out_ <- limma(convertedData.out = convertedData.out_, 
                    readSampleInfo.out = readSampleInfo.out_, 
                    input_dataFileFormat = 1, 
                    input_countsLogStart = 4, 
                    convertedCounts.out = convertedCounts.out_, 
                    input_CountsDEGMethod = 3, 
                    input_limmaPval = 0.1, 
                    input_limmaFC = 2, 
                    input_selectModelComprions = "Gravity: Microgravity vs. Terrestrial", 
                    input_selectFactorsModel = "Gravity", 
                    input_selectInteractions = NULL, 
                    input_selectBlockFactorsModel = NULL, 
                    factorReferenceLevels.out = "Gravity:Terrestrial")
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
DEG.data.out_ <- DEG.data(limma.out = limma.out_, 
                          convertedData.out = convertedData.out_, 
                          allGeneInfo.out = allGeneInfo.out_)

limma.out_$comparisons
## [1] "Microgravity-Terrestrial"
vennPlot(limma.out = limma.out_, 
         input_selectComparisonsVenn = limma.out_$comparisons[1:min(3, length(limma.out_$comparisons))], # if less than three comparisons, include all comparisons.
         input_UpDownRegulated = FALSE) # Split up and down regulated genes

sigGeneStats(limma.out_) # number of DEGs as figure

sigGeneStatsTable(limma.out_) # number of DEGs as table
##                                       Comparisons  Up Down
## Microgravity-Terrestrial Microgravity-Terrestrial 353 1071

7. DEG2

# input_selectContrast_ <- "Microgravity-Terrestrial" # Selected comparisons
selectedHeatmap.data.out_ <- selectedHeatmap.data(convertedData.out = convertedData.out_, 
                                                  readSampleInfo.out = readSampleInfo.out_, 
                                                  limma.out = limma.out_, 
                                                  .converted.out = NULL, 
                                                  .readData.out = readData.out_, 
                                                  .input_noIDConversion = TRUE, 
                                                  input_dataFileFormat = 1, 
                                                  input_CountsDEGMethod = 3, 
                                                  input_selectModelComprions = "Gravity: Microgravity vs. Terrestrial", 
                                                  input_selectFactorsModel = "Gravity", 
                                                  factorReferenceLevels.out = c("Gravity:Terrestrial"), 
                                                  input_selectContrast = "Microgravity-Terrestrial")

selectedHeatmap(selectedHeatmap.data.out = selectedHeatmap.data.out_, 
                .mycolors = mycolors_, 
                .heatColors = heatColors_) # heatmap for DEGs in selected comparison

# Save gene lists and data into files
write.csv(
  selectedHeatmap.data(convertedData.out = convertedData.out_, 
                       readSampleInfo.out = readSampleInfo.out_,
                       .converted.out = converted.out_, 
                       .readData.out = readData.out_, 
                       .input_noIDConversion = input_noIDConversion_, 
                       input_dataFileFormat = 1, 
                       input_CountsDEGMethod = 3, 
                       input_selectModelComprions = "Gravity: Microgravity vs. Terrestrial", 
                       input_selectFactorsModel = "Gravity", 
                       factorReferenceLevels.out = "Gravity:Terrestrial",
                       input_selectContrast = "Microgravity-Terrestrial", 
                       limma.out = limma.out_)$genes, 
  "heatmap.data.csv"
)

write.csv(DEG.data(limma.out = limma.out_, convertedData.out = convertedData.out_, allGeneInfo.out = allGeneInfo.out_), 
          "DEG.data.csv")

write(AllGeneListsGMT(limma.out_), 
      "AllGeneListsGMT.gmt")
input_selectGO2_ <- "GOBP" # Gene set category
geneListData.out_ <- geneListData(limma.out = limma.out_, 
                                  allGeneInfo.out = allGeneInfo.out_, 
                                  input_selectGO2 = "GOBP", 
                                  input_selectOrg = "NEW", 
                                  input_limmaPval = 0.1, 
                                  input_limmaFC = 2, 
                                  input_selectContrast = "Microgravity-Terrestrial")

volcanoPlot(limma.out = limma.out_, 
            input_limmaPval = 0.1, 
            input_limmaFC = 2, 
            input_selectContrast = "Microgravity-Terrestrial")

scatterPlot(limma.out = limma.out_, 
            convertedData.out = convertedData.out_, 
            readSampleInfo.out = readSampleInfo.out_, 
            input_dataFileFormat = 1,
            input_CountsDEGMethod = 3, 
            input_limmaPval = 0.1, 
            input_limmaFC = 2, 
            input_selectModelComprions = "Gravity: Microgravity vs. Terrestrial", 
            input_selectFactorsModel = "Gravity", 
            factorReferenceLevels.out = "Gravity:Terrestrial", 
            input_selectContrast = "Microgravity-Terrestrial")

MAplot(
  limma.out = limma.out_, 
  convertedData.out = convertedData.out_, 
  readSampleInfo.out = readSampleInfo.out_, 
  .converted.out = converted.out_, 
  .readData.out = readData.out_, 
  .input_noIDConversion = TRUE, 
  input_dataFileFormat = 1, 
  input_CountsDEGMethod = 3, 
  input_limmaPval = 0.1, 
  input_limmaFC = 2, 
  input_selectModelComprions = "Gravity: Microgravity vs. Terrestrial", 
  input_selectFactorsModel = "Gravity", 
  factorReferenceLevels.out = "Gravity:Terrestrial", 
  input_selectContrast = "Microgravity-Terrestrial"
)

geneListGOTable.out_ <- geneListGOTable(GeneSets.out = GeneSets.out_, 
                                        minGenesEnrichment = 2, 
                                        selectedHeatmap.data.out = selectedHeatmap.data.out_)
# Read pathway data again
GeneSets.out_ <- readGeneSets(
  fileName = geneSetFile_,
  convertedData = convertedData.out_, 
  GO = "GOBP", 
  selectOrg = "NEW",
  myrange = c(15, 2000)
)

#input_removeRedudantSets_ <- TRUE # Remove highly redundant gene sets?
results <- geneListGO(geneListGOTable.out = geneListGOTable.out_, 
                      input_removeRedudantSets = TRUE) # Enrichment analysis

results$adj.Pval <- format(results$adj.Pval, digits = 3)

kable(results, row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
Direction adj.Pval nGenes Pathways
Down regulated 4.1e-37 143 Response to external stimulus
4.6e-36 117 Response to biotic stimulus
7.5e-35 114 Response to external biotic stimulus
7.5e-35 114 Response to other organism
1.1e-34 124 Defense response
2.3e-33 134 Multi-organism process
1.7e-28 60 Secondary metabolic process
1.1e-27 79 Response to drug
1.5e-26 137 Response to oxygen-containing compound
2.8e-26 156 Response to organic substance
Up regulated 1.9e-18 29 Photosynthesis
5.7e-15 61 Organonitrogen compound biosynthetic process
7.5e-15 41 Peptide metabolic process
5.2e-14 38 Translation
5.2e-14 38 Peptide biosynthetic process
1.6e-13 42 Cellular amide metabolic process
1.6e-13 39 Amide biosynthetic process
3.1e-13 19 Pigment biosynthetic process
1.1e-12 20 Pigment metabolic process
1.3e-11 23 Plastid organization

STRING-db API access. We need to find the taxonomy id of your species, this used by STRING. First we try to guess the ID based on iDEP’s database. Users can also skip this step and assign NCBI taxonomy id directly by findTaxonomyID.out_ = 10090 # mouse 10090, human 9606 etc.

STRING10_species_ <- read.csv(STRING10_speciesFile_)
ix <- grep("Arabidopsis thaliana", STRING10_species_$official_name)
findTaxonomyID.out_ <- STRING10_species_[ix, 1] # find taxonomyID
findTaxonomyID.out_
## [1] 3702

Enrichment analysis using STRING

STRINGdb_geneList.out_ <- STRINGdb_geneList(geneListData.out = geneListData.out_, 
                                            findTaxonomyID.out = findTaxonomyID.out_) # convert gene lists
## Warning:  we couldn't map to STRING 0% of your identifiers
# input_STRINGdbGO_ <- "Process" #' Process', 'Component', 'Function', 'KEGG', 'Pfam', 'InterPro'
results <- stringDB_GO_enrichmentData(selectedHeatmap.data.out = selectedHeatmap.data.out_, 
                                      minGenesEnrichment = 2, 
                                      input_STRINGdbGO = "Process", #' Process', 'Component', 'Function', 'KEGG', 'Pfam', 'InterPro'
                                      findTaxonomyID.out = findTaxonomyID.out_, 
                                      STRINGdb_geneList.out = STRINGdb_geneList.out_) # enrichment using STRING
## Warning in string_db$get_enrichment(ids, category = input_STRINGdbGO, methodMT =
## "fdr", : methodMT parameter is depecated. Only FDR correction is available.
## Warning in string_db$get_enrichment(ids, category = input_STRINGdbGO, methodMT =
## "fdr", : iea parameter is deprecated.
## [1] "Process"
## Warning in string_db$get_enrichment(ids, category = input_STRINGdbGO, methodMT =
## "fdr", : methodMT parameter is depecated. Only FDR correction is available.

## Warning in string_db$get_enrichment(ids, category = input_STRINGdbGO, methodMT =
## "fdr", : iea parameter is deprecated.
## [1] "Process"
results$adj.Pval <- format(results$adj.Pval, digits = 3)

kable(results, row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
“No significant enrichment found.” adj.Pval
No significant enrichment found. NULL

PPI network retrieval and analysis

stringDB_network1(geneLists = 1, 
                  input_nGenesPPI = 100, 
                  findTaxonomyID.out = findTaxonomyID.out_, 
                  STRINGdb_geneList.out = STRINGdb_geneList.out_) # Show PPI network

Generating interactive PPI

write(
  stringDB_network_link(
    input_nGenesPPI = 100, 
    findTaxonomyID.out = findTaxonomyID.out_, 
    STRINGdb_geneList.out = STRINGdb_geneList.out_, 
    geneListData.out = geneListData.out_
  ), 
  "PPI_results.html"
)
## Warning: 'string_db$get_link' is deprecated.
## Use 'Contact developers to request functionality' instead.
## See help("Deprecated")
## Warning:  we couldn't map to STRING 0% of your identifiers
## Warning: 'string_db$get_link' is deprecated.
## Use 'Contact developers to request functionality' instead.
## See help("Deprecated")

## Warning: 'string_db$get_link' is deprecated.
## Use 'Contact developers to request functionality' instead.
## See help("Deprecated")

8. Pathway analysis

# Read pathway data again
GeneSets.out_ <- readGeneSets(
  fileName = geneSetFile_,
  convertedData = convertedData.out_, 
  GO = "GOBP", 
  selectOrg = "NEW",
  myrange = c(15, 2000)
)

gagePathwayData.out_ <- gagePathwayData(limma.out = limma.out_, 
                                        input_minSetSize = 15, 
                                        input_maxSetSize = 2000, 
                                        input_selectContrast1 = "Microgravity-Terrestrial", 
                                        input_pathwayPvalCutoff = 0.2, 
                                        input_nPathwayShow = 30, 
                                        input_absoluteFold = FALSE, 
                                        input_GenePvalCutoff = 1, 
                                        GeneSets.out = GeneSets.out_) # pathway analysis using GAGE

results <- gagePathwayData.out_ # Enrichment analysis for k-Means clusters
results$adj.Pval <- format(results$adj.Pval, digits = 3)
kable(results, row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
Direction GAGE analysis: Microgravity vs Terrestrial statistic Genes adj.Pval
Down Response to chitin -7.7014 104 1.9e-10
Secondary metabolic process -7.4839 265 1.9e-10
Response to drug -7.3469 481 1.9e-10
Response to organonitrogen compound -6.3515 212 1.3e-07
Secondary metabolite biosynthetic process -6.3149 118 2.4e-07
Response to nitrogen compound -6.2461 271 1.6e-07
Up Photosynthesis 10.7757 223 3.7e-21
Ribosome biogenesis 10.2165 342 6.5e-20
Ribonucleoprotein complex biogenesis 10.0934 437 6.5e-20
NcRNA metabolic process 9.4896 425 8.0e-18
Plastid organization 9.039 257 6.2e-16
NcRNA processing 8.7867 357 2.8e-15
RRNA processing 8.2838 239 2.6e-13
RRNA metabolic process 8.2747 244 2.6e-13
Photosynthesis, light reaction 7.9994 119 7.6e-12
Chloroplast organization 7.4927 198 4.2e-11
Generation of precursor metabolites and energy 7.2638 391 7.9e-11
RNA modification 7.1316 321 2.8e-10
Thylakoid membrane organization 6.8664 46 1.5e-07
Tetrapyrrole biosynthetic process 6.0669 70 7.8e-07
Ribosome assembly 6.0176 76 1.3e-06
Ribosomal large subunit biogenesis 5.9801 99 1.0e-06
Tetrapyrrole metabolic process 5.964 93 7.8e-07
Porphyrin-containing compound biosynthetic process 5.9033 67 1.3e-06
Porphyrin-containing compound metabolic process 5.8821 92 1.1e-06
Protein transmembrane transport 5.7744 112 1.4e-06
Ribonucleoprotein complex subunit organization 5.7385 181 1.3e-06
Nucleoside monophosphate metabolic process 5.6492 193 1.6e-06
Ribonucleoprotein complex assembly 5.6323 173 1.8e-06
TRNA metabolic process 5.6148 170 1.8e-06 
pathwayListData.out_ <- pathwayListData(allGeneInfo.out = allGeneInfo.out_, 
                                        input_selectOrg = "NEW", 
                                        input_selectGO = "GOBP", 
                                        input_pathwayMethod = 1, 
                                        gagePathwayData.out = gagePathwayData.out_, 
                                        fgseaPathwayData.out = fgseaPathwayData.out_, 
                                        GeneSets.out = GeneSets.out_)

enrichmentPlot(pathwayListData.out_, enrichedTerms = 25)
## NULL
enrichmentNetwork(pathwayListData.out_)

enrichmentNetworkPlotly(pathwayListData.out_)

# input_pathwayMethod_ <- 3 # 1  fgsea

fgseaPathwayData.out_ <- fgseaPathwayData(limma.out = limma.out_, 
                                          input_minSetSize = 15, 
                                          input_maxSetSize = 2000, 
                                          input_selectContrast1 = "Microgravity-Terrestrial", 
                                          input_pathwayPvalCutoff = 0.2, 
                                          input_absoluteFold = FALSE, 
                                          input_nPathwayShow = 30, 
                                          input_GenePvalCutoff = 1, 
                                          GeneSets.out = GeneSets.out_) # Pathway analysis using fgsea
## Warning in fgsea(pathways = gmt, stats = fold, minSize = input_minSetSize, :
## You are trying to run fgseaSimple. It is recommended to use fgseaMultilevel. To
## run fgseaMultilevel, you need to remove the nperm argument in the fgsea function
## call.
## Warning in fgseaSimple(...): There were 4 pathways for which P-values were not
## calculated properly due to unbalanced gene-level statistic values
results <- fgseaPathwayData.out_ # Enrichment analysis for k-Means clusters
results$adj.Pval <- format(results$adj.Pval, digits = 3)
kable(results, row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
Direction GSEA analysis: Microgravity vs Terrestrial NES Genes adj.Pval
Up Photosynthesis 3.3134 223 1.2e-02
Photosynthesis, light reaction 3.0444 119 6.2e-03
Thylakoid membrane organization 2.9683 46 4.0e-03
Translation 2.9257 619 1.5e-01
Plastid organization 2.9112 257 1.5e-02
Peptide biosynthetic process 2.8963 622 1.5e-01
Photosynthetic electron transport chain 2.854 46 4.0e-03
Ribosome biogenesis 2.8514 342 2.5e-02
Chloroplast organization 2.7515 198 1.0e-02
Ribosome assembly 2.7444 76 4.7e-03
Plastid membrane organization 2.7133 49 4.1e-03
Chlorophyll biosynthetic process 2.7088 58 4.2e-03
Protein localization to chloroplast 2.7023 45 4.0e-03
Tetrapyrrole biosynthetic process 2.6956 70 4.6e-03
RRNA processing 2.6864 239 1.4e-02
Ribonucleoprotein complex biogenesis 2.6828 437 4.0e-02
RRNA metabolic process 2.679 244 1.4e-02
Ribosomal large subunit biogenesis 2.6547 99 5.3e-03
Porphyrin-containing compound biosynthetic process 2.6546 67 4.4e-03
Protein targeting to chloroplast 2.6538 43 3.9e-03
Establishment of protein localization to chloroplast 2.6538 43 3.9e-03
Tetrapyrrole metabolic process 2.6537 93 5.1e-03
Porphyrin-containing compound metabolic process 2.6332 92 5.1e-03
Chlorophyll metabolic process 2.6244 81 4.8e-03
Chloroplast RNA processing 2.5756 19 3.6e-03
Pigment biosynthetic process 2.5662 129 6.6e-03
NcRNA processing 2.5179 357 2.7e-02
Chloroplast rRNA processing 2.5112 18 3.6e-03
Protein transmembrane transport 2.5082 112 5.8e-03
Starch metabolic process 2.5005 60 4.3e-03 
pathwayListData.out_ <- pathwayListData(allGeneInfo.out = allGeneInfo.out_, 
                                        input_selectOrg = "NEW", 
                                        input_selectGO = "GOBP", 
                                        input_pathwayMethod = 3, 
                                        gagePathwayData.out = gagePathwayData.out_, 
                                        fgseaPathwayData.out = fgseaPathwayData.out_, 
                                        GeneSets.out = GeneSets.out_)

enrichmentPlot(enrichedTerms = pathwayListData.out_, 
               rightMargin = 25, 
               mycolors = mycolors_)

enrichmentNetwork(pathwayListData.out_)

enrichmentNetworkPlotly(pathwayListData.out_)

PGSEAplot(convertedData.out = convertedData.out_, 
          readSampleInfo.out = readSampleInfo.out_, 
          input_selectOrg = "NEW", 
          input_dataFileFormat = 1, 
          input_selectGO = "GOBP", 
          input_minSetSize = 15, 
          input_maxSetSize = 2000, 
          input_CountsDEGMethod = 3, 
          input_selectModelComprions = "Gravity: Microgravity vs. Terrestrial", 
          input_selectFactorsModel = "Gravity", 
          factorReferenceLevels.out = "Gravity:Terrestrial", 
          input_selectContrast1 = "Gravity:Terrestrial", 
          input_pathwayPvalCutoff = 0.2, 
          input_nPathwayShow = 30, 
          GeneSets.out = GeneSets.out_) # pathway analysis using PGSEA
## 
## Computing P values using ANOVA

9. Chromosome

#input_selectContrast2_ <- "Microgravity-Terrestrial" # select Comparison
# input_selectContrast2 = limma.out_$comparisons[3] # manually set
#input_limmaPvalViz_ <- 0.1 # FDR to filter genes
#input_limmaFCViz_ <- 2 # FDR to filter genes
genomePlotly(limma.out = limma.out_, 
             allGeneInfo.out = allGeneInfo.out_, 
             input_selectContrast2 = "Microgravity-Terrestrial", 
             input_limmaPvalViz = 0.1, 
             input_limmaFCViz = 2) # shows fold-changes on the genome
## Warning in eval(quote(list(...)), env): NAs introduced by coercion
## Warning in genomePlotly(limma.out = limma.out_, allGeneInfo.out =
## allGeneInfo.out_, : NAs introduced by coercion

10. Biclustering

biclustering.out_ <- biclustering(convertedData.out = convertedData.out_, 
                                  input_nGenesBiclust = 1000, 
                                  input_biclustMethod = "BCCC()") # run analysis

input_selectBicluster_ <- 1 # select a cluster
biclustHeatmap(biclustering.out = biclustering.out_, 
               heatColors = heatColors_, 
               input_heatColors1 = 1, 
               input_selectBicluster = 1) # heatmap for selected cluster

#input_selectGO4_ <- "GOBP" # Gene set category
# Read pathway data again
GeneSets.out_ <- readGeneSets(
  fileName = geneSetFile_,
  convertedData = convertedData.out_, 
  GO = "GOBP", 
  selectOrg = "NEW",
  myrange = c(15, 2000)
)

results <- geneListBclustGO(minGenesEnrichment = 2, 
                            input_selectBicluster = 1, 
                            biclustering.out = biclustering.out_, 
                            GeneSets.out = GeneSets.out_) # Enrichment analysis for k-Means clusters

results$adj.Pval <- format(results$adj.Pval, digits = 3)

kable(results, row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
adj.Pval Genes Pathways
2.4e-134 308 Response to abiotic stimulus
5.9e-86 170 Response to inorganic substance
1.1e-66 216 Response to organic substance
5.7e-64 182 Oxidation-reduction process
1.3e-63 197 Organonitrogen compound biosynthetic process
2.2e-63 173 Response to external stimulus
4.5e-62 144 Response to external biotic stimulus
4.5e-62 144 Response to other organism
2.6e-61 144 Response to biotic stimulus
3.1e-61 105 Response to metal ion

11. Co-expression network

wgcna.out_ <- wgcna(convertedData.out = convertedData.out_, 
                    maxGeneWGCNA = 2000, 
                    input_mySoftPower = 5, 
                    input_nGenesNetwork = 1000, 
                    input_minModuleSize = 20) # run WGCNA
##    Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
## 1      1   0.6880  3.020          0.962  418.00    426.00  557.0
## 2      2   0.5030  1.310          0.949  240.00    241.00  388.0
## 3      3   0.2260  0.549          0.896  159.00    155.00  294.0
## 4      4   0.0337  0.163          0.868  113.00    108.00  233.0
## 5      5   0.0639 -0.187          0.865   85.00     79.00  191.0
## 6      6   0.2370 -0.350          0.875   66.30     60.10  159.0
## 7      7   0.4190 -0.473          0.908   53.20     46.40  135.0
## 8      8   0.5660 -0.576          0.946   43.70     37.00  117.0
## 9      9   0.6550 -0.661          0.944   36.50     30.10  103.0
## 10    10   0.7000 -0.706          0.948   30.90     24.80   91.1
## 11    12   0.7590 -0.822          0.928   23.00     17.50   73.4
## 12    14   0.7970 -0.918          0.922   17.70     12.80   60.7
## 13    16   0.8360 -0.968          0.932   14.00      9.68   51.3
## 14    18   0.8470 -1.020          0.929   11.40      7.55   44.1
## 15    20   0.8340 -1.050          0.905    9.39      5.96   38.5
## TOM calculation: adjacency..
## ..will not use multithreading.
##  Fraction of slow calculations: 0.000000
## ..connectivity..
## ..matrix multiplication (system BLAS)..
## ..normalization..
## ..done.
softPower(wgcna.out_) # soft power curve

listWGCNA.Modules.out <- listWGCNA.Modules(wgcna.out_) # modules
modulePlot(wgcna.out_) # plot modules

# Read pathway data again
GeneSets.out_ <- readGeneSets(
  fileName = geneSetFile_,
  convertedData = convertedData.out_, 
  GO = "GOBP", 
  selectOrg = "NEW",
  myrange = c(15, 2000)
)

moduleNetwork(wgcna.out = wgcna.out_, 
              input_noIDConversion = TRUE, 
              allGeneInfo.out = allGeneInfo.out_, 
              input_selectOrg = "NEW", 
              input_selectWGCNA.Module = "Entire network", 
              input_topGenesNetwork = 10, 
              input_edgeThreshold = 0.4, 
              input_selectGO5 = "GOBP") # show network of top genes in selected module
##  softConnectivity: FYI: connecitivty of genes with less than 4 valid samples will be returned as NA.
##  ..calculating connectivities.. ..0% ..100%

results <- networkModuleGO(GeneSets.out = GeneSets.out_, 
                           minGenesEnrichment = 2, 
                           input_selectWGCNA.Module = "Entire network", 
                           wgcna.out = wgcna.out_) # Enrichment analysis of selected module

results$adj.Pval <- format(results$adj.Pval, digits = 3)
kable(results, row.names = FALSE) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  scroll_box(width = "100%")
adj.Pval Genes Pathways
2.4e-134 308 Response to abiotic stimulus
5.9e-86 170 Response to inorganic substance
1.1e-66 216 Response to organic substance
5.7e-64 182 Oxidation-reduction process
1.3e-63 197 Organonitrogen compound biosynthetic process
2.2e-63 173 Response to external stimulus
4.5e-62 144 Response to external biotic stimulus
4.5e-62 144 Response to other organism
2.6e-61 144 Response to biotic stimulus
3.1e-61 105 Response to metal ion